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ABSTRACT  The  numerical  algorithms  developed  under  the 
6.1  Structural  Modeling  Program  are  reviewed  and  evaluated 
I  to  assess  their  applicability  towards  the  6.2  Nonlinear 

Structural  Analysis  Program. 

j  An  exposition  of  the  currently  available  numerical  meth- 

i  ods  for  the  solution  of  soil-structure  interaction  problems  is 

also  performed.  Within  this  framework  the  work  done  under 
j  the  6.1  Program  is  evaluated. 
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1.  Introduction 

‘  The  6.1  Structural  Modeling  Program  at  the  Naval  Civil  Engineering  Laboratory  at 
Port  Hueneme  has  as  a  main  objective  to  develop  analysis  algorithms  for  the 
nonlinear  analysis  of  drydocks  and  other  large  geotechnical-structural  systems.  The 
6.2  Nonlinear  Structural  Analysis  Program  is  armed  at  implementing  those  algorithms 
into  a  special  purpose  finite  element  analysis  program  that  will  allow  a  3  dimensional 
nonlinear  dynamic  analysis  of  the  large  structural  systems  under  earthquake  and  blast 
type  of  excitations.  This  report  contains  a  review  and  evaluation  of  these  algorithms  to 
assess  their  applicability  towards  the  6.2  Program. 

The  reviewed  work  is  contained  in  References  [1-8].  The  issues  covered  include: 
reduced  modal  methods  for  nonlinear  dynamic  problems,  contact  and  friction  finite 
element  techniques,  combination  of  finite  element  and  recursive  boundary  element 
methods,  validation  of  bounding  surface  plasticity  models  for  saturated  cohesive  soils 
and,  finally,  linear  fracture  mechanics  for  concrete. 

'work  is  initiated  with  an  exposition  of  the  currently  available 
numerical  methods  for  the  solution  of  the  soil-structure  interaction  problems  in  the  time 
and  frequency  domains.  Directions  are  given  as  to  which  of  the  available  methods 
seem  to  be  more  applicable  for  the  problem  at  hand.  In  particular,  complete  and 
substructure  methods  in  the  time  domain  are  presented  as  the  most  convenient  ones. 
In  view  of  these  algorithms  the  work  done  under  the  6.1  Program  is  reviewed  and 
evaluated. 


2.  Numerical  methods  for  soil-structure  interaction  problems 

Although  not  a  part  of  the  6.1  Program  it  is  important  to  explain  and  evaluate  the 
numerical  environment  in  which  the  soil-structure  interaction  problem  involved  in  the 
analysis  of  drydocks  under  earthquake  excitation  and  blast  effects,  is  to  be  analyzed.  A 
number  of  different  methods,  implemented  in  currently  available  computer  codes,  have 
been  used  so  far  for  the  analysis  of  such  problems.  All  these  methods  can  be 
classified  in  two  major  groups:  frequency  domain  and  time  domain  methods. 
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The  first  group  encompasses  the  majority  of  available  methods  and  computer 
codes.  The  available  techniques  within  the  frequency  domain  are  summarized  in 
Fig.1.  These  can  be  classified  into  three  major  groups:  complete,  substructure  and 
hybrid  methods.  Four  different  problems  need  be  solved  for  a  complete  analysis:  the 
site  response  problem  or  free  field  analysis,  the  scattering  of  waves  due  to  the 
presence  of  the  foundation,  the  impedance  problem  that  yields  the  complex  stiffness 
terms  and  the  final  structural  analysis.  The  flexible  volume  methods  seem  to  be  the 
more  attractive  of  all  of  these  because  they  eliminate  the  scattering  problem  and  also 
because  they  take  full  advantage  of  the  substructure  method.  Details  of  the  equations 
involved  in  each  technique  are  given  in  [9,10].  The  common  denominator  in  all  these 
methods  is  the  basic  assumption  of  linearity  in  the  total  response  of  the  system. 

The  reasons  for  the  popularity  of  the  frequency  domain  approaches  are,  firstly, 
the  possibility  of  dividing  the  problem  into  substructures  that  can  be  analyzed 
independently,  and  secondly,  the  availability  of  accurate  frequency  dependent 
radiation  boundaries  that  help  to  considerably  reduce  the  size  of  the  finite  element 
models.  However,  as  pointed  out  in  [9],  these  two  advantages  can  not  hinder  the 
tremendous  limitations  that  arise  from  the  fact  that  frequency  domain  techniques  can 
not  solve  nonlinear  (geometric  or  material)  problems.  Quasi  nonlinear  approaches 
[1 1 ,12]  have  been  developed  that  require  a  continuous  transfer  of  information  between 
the  frequency  and  time  domains.  These  quasi  nonlinear  methods  lack  rigor  and 
demand  an  exceedingly  large  computational  effort.  Another  important  limitation  comes 
from  the  fact  that  frequency  domain  techniques  become  numerically  inefficient  for 
transient  analysis  of  three  dimensional  problems.  In  view  of  these  facts,  frequency 
domain  methods  do  not  seem  the  most  appropriate  for  the  kind  of  problems  involved  in 
the  transient  dynamic  analysis  of  drydocks  that  involve  nonlinear  material  and  contact 
problems. 

A  classification  analogous  to  that  done  with  the  frequency  domain  can  be  done 
with  the  time  domain  methods.  Fig.  2  illustrates  the  available  methods  in  the  time 
domain.  Three  major  groups  are  distinguished:  the  complete,  the  boundary  and  the 
volume  methods.  There  is  no  such  a  thing  as  an  impedance  problem,  however  the 
structural  analysis  now  involves  the  total  ensemble  soil-structure.  In  this  regard  it  is 
important  to  note  that  the  substructure  concept  used  in  the  frequency  domain  is  not 
similar  to  that  currently  used  in  the  time  domain.  The  splitting  of  the  model,  as  done  in 
the  frequency  domain,  is  substituted  in  the  time  domain  by  a  reduction  in  the  number 
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of  degrees  of  freedom  in  parts  of  the  system  (substructures)  that  are  subsequently 
assembled  and  solved  simultaneously.  An  important  observation  is  that  the 
mechanisms  available  in  the  time  domain  to  account  for  the  radiation  of  energy 
through  the  boundaries  of  the  mesh  are  not  as  accurate  as  those  available  in  the 
frequency  domain.  This  limitation  is  compensated  by  having  recourse  to  larger  finite 
element  models,  and  approximate  frequency  independent  dashpots.  Nevertheless, 
satisfactory  solutions  to  the  radiation  problem  in  the  time  domain  have  been  obtained 
as  reported  in  References  [9,13,14], 

The  volume  methods  are  not  applicable  to  the  drydock  problem  since  it  involves  a 
superposition  scheme  at  the  foundation  level  which  obviously  precludes  any  nonlinear 
analysis.  The  two  other  possibilities  are  the  complete  method  and  the  boundary 
method.  The  first  one  allows  a  full  nonlinear  analysis  of  all  the  ensemble  soil-structure 
problem,  however  it  requires  an  extensive  modeling  of  the  soil  profile  until  a 
sufficiently  rigid  bedrock,  at  which  the  input  motion  is  to  be  specified,  is  encountered. 
This  suffers  from  the  major  drawback  that  unless  the  bedrock  is  close  to  the  surface  the 
3-dimensional  model  needed  for  the  analysis  may  become  extremely  large,  thus 
making  the  computational  cost  of  analysis  prohibitive. 

In  the  boundary  method  the  input  motion  can  be  specified  at  a  surface  marked  by 
the  separation  between  the  linear  and  nonlinear  soil  characteristics,  as  shown  in  Fig. 
3.  Its  major  advantages  lie  in  the  following  facts: 

-  The  soil  below  the  input  boundary  can  be  considered  linear  and  therefore  that 
part  of  the  finite  element  model  may  be  substructured,  thus  reducing  the  size  of 
the  final  model.  This  same  substructure  can  be  used  to  solve  the  scattering 
problem  avoiding  part  of  the  extra  computational  effort  involved  in  this  phase  of 
the  problem. 

-  The  soil  profile  can  be  cut  at  a  closer  distance  to  the  structure  than  that  of  the 
complete  methods.  This  distance  will  exclusively  depend  on  the  accuracy 
required  from  the  energy  radiation  mechanism  at  the  boundaries  of  the  finite 
element  model. 

On  the  other  hand  the  boundary  method  has  the  disadvantage,  as  compared  to 
the  complete  method,  of  having  to  solve  the  scattering  problem.  The  input  motion, 
although  of  a  multiple  excitation  type,  is  applied  in  a  smaller  amount  of  nodes  than  in 
the  case  of  the  complete  methods. 


The  factors  mentioned  above  should  be  weighted  with  the  typical  properties  of 
the  drydock  environment  and  soil  characteristics  to  decide  which  of  the  two  proposed 
methods  should  be  selected  for  the  special  purpose  analysis  program  under 
consideration. 

There  is  no  doubt  that  a  3-D  analysis  of  the  types  of  drydocks  under  consideration 
will  involve  a  much  greater  computational  effort  than  a  2-D  analysis.  However,  it  is 
important  to  point  out  that  the  former  becomes  necessary,  first,  when  the  loading  has 
strong  components  in  the  normal  direction  and,  second,  when  the  structural  system 
contains  parts,  sometimes  of  vital  importance,  which  do  not  have  plane  symmetry  and 
therefore  need  be  analyzed  under  3  dimensional  assumptions.  In  particular,  for  the 
drydock  problem,  the  analysis  of  the  pump  house  and  the  caisson  require  full  3-D 
consideration. 

With  these  ideas  in  mind  we  can  now  review  the  work  done  under  the  6.1 
Program  and  see  how  they  fit  into  the  general  formulation  of  the  problem. 

3.  Reviewed  Methods 

In  the  present  study  several  methods  developed  under  the  6.1  Structural 
Modeling  Program  have  been  reviewed  to  determine  their  applicability  towards  the 
objectives  of  the  6.2  Nonlinear  Structural  Analysis  Program  aimed  at  drydock  analysis. 
The  methods  considered  include: 

a)  Reduced  coordinates  for  the  solution  of  large  time  dependent  problems. 

b)  Contact  problems. 

c)  Combined  boundary  element  and  finite  element  method  for  the  silent 
boundary  effect. 

d)  Bounding  surface  plasticity  model  for  soils,  and  linear  elastic  fracture 
mechanics  for  concrete. 

3.1  REDUCED  COORDINATES  FOR  LARGE  TIME  DEPENDENT  PROBLEMS 

The  work  reviewed  in  this  section  is  reported  in  References  [1,2].  This  work  is 
mainly  devoted  to  the  construction  of  algorithms  that  will  reduce  the  size  of  finite 


element  equations  and  will  retain  the  fundamental  properties  of  the  larger  original 
system. 

The  proposed  method  can  be  basically  summarized  as  follows: 

1.  Use  of  a  series  of  "Lanczos  vectors"  (generated  taking  into  account  the  spatial 
distribution  of  the  loading)  or  "Lanczos  eigenvectors"  to  reduce  the  total  size  of 
the  finite  element  equations.  The  vectors  are  originated  in  an  inverse  setting  so 
that  the  projection  basis  is  related  to  the  smallest  eigenpairs  of  the  system.  If 
resonance  is  expected,  shifting  may  be  used  to  generate  vectors  in  a  given 
interval  of  the  frequency  spectrum. 

2.  Use  a  Newmark-Newton  time  marching  scheme  for  the  integration  of  the 
reduced  set  of  equations.  Each  step  will  require  the  transformation  of  the 
nonlinear  stiffness  matrices  from  the  finite  element  coordinates  to  the  Lanczos 
coordinates,  operations  that  can  be  performed  at  the  element  level,  and  the 
solution  of  the  reduced  system  of  equations.  Two  schemes  are  proposed  to 
assure  convergence.  The  first  one  is  a  time  sub-stepping  with  an  integer  number 
of  substeps.  The  second  is  a  recalculation  of  the  vector  basis,  which  will  be 
necessary  if  the  finite  element  mesh  undergoes  large  deformations.  The 
updating  of  the  basis  constitutes  one  of  the  most  computationally  expensive 
parts  of  the  procedure  and  this  cost  may  be  critical  for  the  overall  efficiency  of  the 
method:  If  neither  of  these  two  schemes  allows  convergence,  the  algorithm  will 
have  to  be  stopped.  Alternatives  to  this  problem  are  currently  being  searched  for 
by  the  authors. 

The  method  has  been  tested  in  several  linear  and  nonlinear  problems  yielding 
very  good  results. 

The  proposed  method  should  be  very  applicable  to  the  nonlinear  analysis  of  the 
drydock.  No  testing  of  the  algorithm  in  large  scale  problems  have  been  reported, 
however  all  the  indicated  operation  counts  point  to  the  fact  that  this  reduction  method 
should  overcome  the  direct  methods  in  computational  efficiency.  This  algorithm  will 
also  allow  the  analysis  of  the  large  drydock  finite  element  model  in  a  non¬ 
supercomputer  environment.  Several  considerations  may  be  suggested: 

a)  The  use  of  substructures  for  the  analysis  at  hand  is  almost  imperative.  Major 

advantages  in  the  use  of  substructures  may  be  found  when  the  same  structural 


element  is  repeated  several  times,  and  in  the  possibility  of  analysis  of  different 
components  by  different  teams.  In  addition,  as  explained  by  Clough  and  Wilson 
in  [15),  linear  substructures  do  not  need  to  be  updated  during  the  nonlinear 
iteration  process.  The  use  of  substructures  and  mode  superposition  in  nonlinear 
systems  can  also  be  done  as  shown  by  Bathe  and  Gracewski  in  [16].  Each  of 
the  main  physical  components  of  the  problem:  caisson,  pumphouse,  drydock 
walls,  near  field  soil  and  far  field  soil,  can  be  substructured  individually  using 
classical  techniques  [17-19],  The  number  of  boundary  (master)  degrees  of 
freedom  can  be  very  large  (particularly  in  the  walls  of  the  drydock  and  in  the  soil 

t 

system),  in  this  case  a  multilevel  substructuring  technique  (define  substructures 
of  substructures)  may  be  used  as  explained  by  Wilson  and  Bayo  in  [20]. 
Proceeding  in  this  fashion  will  minimize  the  use  of  back-up  storage,  and  will 
allow  for  the  solution  of  the  final  set  of  equations  directly  in  high-speed  storage. 

A  transformation  of  the  component  synthesis  type  can  be  applied  at  each 
substructure  level: 

fi  Tj  1  [Yi 

0  1  Ub  (1) 

where  Uj  and  ub  represent  the  internal  and  boundary  degrees  of  freedom 
respectively,  fj  and  Yj  are  the  Lanczos  vectors  and  generalized  coordinates, 
respectively.  I  is  the  identity  matrix  and  Tj  is  the  static  condensation 
transformation:  -  [Kjj]-1  [Kib]  which  can  be  obtained  by  backsubstitution  of  [Kjj]  [Tj] 
=  -[Kit].  The  triangularization  of  Kjj  can  also  be  used  for  the  calculation  of  the 
Lanczos  vectors.  The  proposed  algorithm  for  the  generation  of  reduced 
Lanczos  coordinates  can  then  be  used  at  each  substructure  level  and  in  the 
global  system  to  reduce  their  size,  and  to  capture  the  inertia  effects  of  the 
internal  degrees  of  freedom.  This  substructuring  part  will  require  serious 
thought  at  the  time  of  computer  implementation. 

b)  As  reported  in  [20,21]  the  use  of  the  "Lanczos  vectors",  which  are  generated 
taking  into  consideration  the  spatial  distribution  of  the  loading,  gives  better  (or 
at  least  as  good)  results  for  earthquake  excitation  than  do  the  "Lanczos 
eigenvectors".  Furthermore  the  former  are  generated  at  a  fraction  of  the  cost 


involved  in  obtaining  the  structural  eigenvectors.  It  is  therefore  recommended 
that  the  "Lanczos  vectors"  be  used  instead  of  the  "Lanczos  eigenvectors". 


c)  A  general  3-D  earthquake  excitation  with  six  different  components  can  be 
expressed  as  a  product  of  spatial  load  vectors  and  time  functions,  as  follows: 


o 


where  fj(s)  represents  the  ith  component  of  the  spatial  distribution  of  the  loading 
and  g,(t)  the  ith  time  variation  function.  The  Krylov  subspace  can  be  generated 
in  this  case  by  using  a  recurrence  sequence  that  uses  "blocks"  of  vectors 
corresponding  initially  to  the  multiload  patterns.  It  is  suggested  that  this  type  of 
generation  be  used  in  the  computer  program.  Methods  to  generate  the  block 
Lanczos  vectors  are  presented  in  [21-23]. 

d)  It  may  be  more  computationally  efficient  to  use  the  selective 
reorthogonalization  scheme  proposed  by  Simon  in  [24]  and  used  by  Nour-Omid 
et  al  in  [25,26],  than  the  one  presented  by  Mish  in  the  report  under 
consideration.  Also,  the  one  proposed  by  Leger  in  [21]  yields  very  accurate 
results  and  is  not  so  computationally  expensive. 

e)  As  a  final  comment  it  is  worth  pointing  out  that  the  proposed  Newton-Lanczos 
algorithm  should  be  a  good  tool  for  the  3-D  nonlinear  analysis  of  drydock  and 
given  the  large  size  of  the  model,  the  only  one  which  can  make  the  problem 
tractable  for  the  computer  environment  available  at  the  NCEL.  However,  the 
programing  effort  involved  in  the  computer  implementation  of  the  proposed 
Lanczos  algorithm  within  the  context  of  nonlinear  analysis  and  multilevel 
substructuring  is  not  an  easy  task. 

3.2  CONTACT  PROBLEMS 

The  work  reviewed  in  this  section  is  reported  in  References  [3-4],  This  work  deals 
with  the  construction  of  algorithms  which  will  handle  dynamic  contact  simulations. 
Landers  in  [3]  reviews  several  frictionless,  small  deformation  contact  algorithms,  and 
discusses  their  advantages  and  disadvantages.  The  conclusion  of  the  study  is  that  the 
augmented  Lagrangian  formulation  combines  the  best  features  of  the  other  currently 
available  methods,  namely,  the  penalty,  the  classical  Lagrangian  and  the  perturbed 
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Lagrangian  methods.  This  formulation  yields  an  algorithm  that  is  not  as  sensitive  to  the 
penalty  number  as  the  penalty  and  perturbed  Lagrangian  methods,  and  does  not 
increase  the  number  of  equations  as  the  classical  Lagrangian  approach. 

In  addition  the  problem  of  the  finite  element  system  of  equations  being  indefinite 
is  also  avoided.  In  the  augmented  Lagrangian  algorithm  the  Lagrange  multiplier  is 
simply  added  to  the  right  hand  side  of  the  equation  as  a  correction  term  that  is 
calculated  through  an  iteration  process.  Convergence  is  assured  without  requiring  that 
the  penalty  parameter  tend  to  infinity.  Also,  the  displacements  and  contact  forces  are 
obtained  without  the  need  for  extra  calculations.  The  work  done  by  Ju  el  al  in  [4]  is 
concerned  with  the  development  of  a  perturbed  Lagrangian  formulation  that  includes 
friction  forces.  The  formulation  is  based  on  an  operator  split  methodology  similar  to 
that  used  in  elasto-plasticity. 

Both  numerical  algorithms  have  been  implemented  in  the  computer  program 
FEAP,  and  tested  through  the  solution  of  different  sample  problems,  yielding  very  good 
results.  Therefore,  they  seem  very  suitable  for  the  analysis  of  the  contact  problem 
existing  between  the  walls  of  the  drydock  and  the  adjacent  soil. 

However,  there  are  at  present  some  serious  limitations  in  the  finite  element 
implementation  of  both  techniques.  In  particular,  the  contact  kinematic  assumptions 
used  in  [3-4]  are  only  currently  valid  for  2-D  problems,  single  contact  segments  and 
linear  interpolation  elements.  The  extension  to  3-D  cases,  multiple  contact  segments, 
symmetric  boundary  conditions  and  higher  order  interpolation  is  still  a  matter  of  basic 
research.  Furthermore,  the  computational  cost  required  to  perform  the  penetration 
calculations  in  the  3-D  problem  might  be  larger  than  that  required  to  solve  the  entire 
finite  element  system  of  equations. 

In  order  to  make  the  problem  tractable  at  the  current  stage,  simpler  kinematic 
assumptions  should  be  made  at  the  time  of  implementing  the  proposed  algorithms  in 
the  computer  program  for  the  analysis  of  the  drydock.  It  is  suggested  that  for  the  first 
version  of  the  program,  small  deformation  and  linear  kinematics  based  on  elementary 
node-on-node  contact  or  a  very  simple  non-nodal  contact  should  be  used.  These 
assumptions  were  used  in  less  recent  work  in  contact  problems,  and  are  reported  in 
[27-28].  These  kinematic  hypotheses  could  be  improved  and  updated,  within  the  finite 
element  code,  as  research  in  the  area  of  3-D  contact  problem  progresses. 
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3.3  COMBINED  BOUNDARY  ELEMENT  AND  FINITE  ELEMENT  METHODS  FOR  THE 
SILENT  BOUNDARY  EFFECT 

The  work  reviewed  in  this  section  is  reported  in  Reference  [5].  The  ultimate 
objective  of  this  research  is  to  combine  the  finite  element  method  with  the  boundary 
element  method  in  an  efficient  manner  for  the  solution  of  large  structural-geotechnical 
problems.  It  is  widely  acknowledged  that  the  former  permits  a  detailed  modeling  of  a 
structural  or  soil  system  and  provides  a  suitable  computational  environment  for  the 
solution  of  dynamic  and  nonlinear  problems.  However,  as  pointed  out  above,  up  to  this 
date  there  are  no  accurate  solutions  within  the  finite  element  framework  to  problems 
involving  the  modeling  of  an  infinite  or  semi-infinite  media  in  the  time  domain.  Besides, 
the  polynomial  expansions  used  in  the  finite  element  method  fail  to  model  singularities 
(at  which  all  or  part  of  the  derivatives  are  infinite)  in  an  accurate  way. 

On  the  other  hand,  the  boundary  element  method  requires  only  the  boundary  to 
be  discretized  and  is  based  on  the  Green's  function  solution  for  the  rest  of  the  domain. 
This  method  permits  a  better  analysis  of  infinite  domain  and  singular  problems.  Due  to 
the  fact  that  this  technique  has  not  been  as  extensively  developed  as  the  finite  element 
method  for  nonlinear-dynamic  problems,  it  seems  that  a  marriage  between  both  could 
provide  appropriate  general  solutions  to  a  large  number  of  problems  in 
elastodynamics,  and  in  particular  to  the  soil-structure  interaction  problem  involved  in 
the  drydock  analysis. 

Previous  studies  at  the  NCEL  have  focussed  on  the  coupling  between  the 
boundary  element  and  the  finite  element  models,  and  also  on  the  improvement  of  the 
performance  of  the  boundary  elements  in  the  near  boundary  region.  The  first  issue  has 
been  addressed  by  Shugar  and  Cox  in  [29],  where  successful  results  have  been 
obtained  in  two-dimensional  elastostatic  problems  with  constant  stress  elements. 
However,  the  fact  that  the  resultant  stiffness  matrix  is  unsymmetric  did  not  make  the 
method  suitable  for  application  in  nonlinear  soil-structure  interaction  problems. The 
reviewed  report  dedicates  attention  to  procedures  that  will  increase  the  reliability  of  the 
method  in  the  near-boundary  region  where  error  is  normally  excessive,  and  in  this 
respect  the  authors  devise  a  recursive  element  subdivision  procedure  for  both  indirect 
and  direct  boundary  element  methods,  that  yields  improved  accuracy  while 
maintaining  simplicity  in  the  implementation.  The  proposed  adaptive  method  can  also 
be  extended  to  the  combination  boundary  and  finite  elements,  but  this  part  needs 
further  research. 
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Unfortunately,  techniques  to  improve  the  matching  between  both  the  finite 
element  and  boundary  element  methods,  as  well  as  the  finding  of  general  integral 
solutions  for  the  three-dimensional  far  field  problem  in  the  time  domain  are  still  under 
basic  research,  and  at  this  stage  they  can  not  be  applied  in  the  6.2  Nonlinear 
Structural  Analysis  Program.  More  specifically,  further  work  is  needed  to  produce  a 
stiffness  coupling  between  the  boundary  element  and  the  finite  element  methods  that, 
for  higher  order  interpolations,  does  not  violate  compatibility  along  the  boundary  and 
that  produces  symmetric  stiffness  matrices.  The  problem  of  finding  a  solution  to  the  3-D 
far  field  problem  could  be  somehow  connected  with  the  problem  of  using  the  boundary 
element  method  to  develop  infinite  elements  in  the  time  domain. 

In  order  to  cope  with  the  problem  of  the  silent  boundary  effect,  it  is  suggested  that 
finite  element  techniques  in  the  time  domain,  reported  above,  be  used  in  conjuction 
with  the  proposed  substructuring  procedure  to  account  for  the  energy  dissipation  at  the 
boundaries  of  the  finite  element  model. 


3.4  BOUNDING  SURFACE  PLASTICITY  MODEL  FOR  SOILS  AND  LINEAR  ELASTIC 
FRACTURE  MECHANICS  FOR  CONCRETE 

References  [6-7]  describe  the  work  done  under  the  6.1  Structural  Modeling 
Program  regarding  the  implementation  and  validation  of  the  bounding  surface 
plasticity  model  for  cohesive  soils.  The  first  work  describes  the  use  of  numerical 
algorithms  to  correct  certain  computational  problems  of  previous  implementations, 
which  already  had  the  capability  to  model  real  soil  behavior.  The  computational 
problems  were  basically  the  numerical  integration  error  occurring  when  using 
moderate  sized  integration  steps,  and  the  inadequacy  of  the  scaling  procedure  to 
return  a  predicted  stress  state  to  the  bounding  surface  when  it  fell  outside. 


Robustness  has  been  added,  firstly  through  the  use  of  sub-stepping  that  allows 
for  convergence  even  in  the  case  of  large  size  integration  steps,  and  secondly, 
through  a  more  accurate  radial  return  algorithm  to  bring  a  point  back  to  the  bounding 
surface.  The  improved  bounding  surface  plasticity  model  is  included  in  the  master 
subroutine  CLAY,  which  can  be  added  to  existing  finite  element  codes  to  supply  the 
incremental  material  properties  needed  in  the  nonlinear  analysis.  The  CLAY 
subroutine  package  allows  the  use  of  a  global  iteration  scheme,  Newton-Raphson  or 
quasi-Newton  methods  for  non-linear  iteration.  It  also  contemplates  the  use  of  a  nearly 
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incompressible  solid  model  for  the  soil  for  those  situations  in  which  the  soil  is 
saturated  in  undrained  conditions. 


The  analysis  program  has  been  validated  as  reported  in  [7]  by  means  of  a 
comparison  with  a  complete  centrifuge  test  of  a  retaining  wall  system.  The  good 
correlations  obtained  between  the  experimental  data  and  the  finite  element  analysis 
corroborate  both  the  reliability  of  the  plasticity  model  and  the  accuracy  of  the  CLAY 
computer  package. 


Although  the  impact  that  the  modifications  and  the  material  model  evaluation 
have  on  the  total  computational  cost  for  large  finite  element  analysis  have  not  yet  been 
determined,  it  may  be  concluded  that  the  technique  is  ready  for  implementation  in  the 
6.2  Nonlinear  Structural  Analysis  Program. 

Reference  [8]  describes  the  work  done  under  the  6.1  Program  to  find  an 
appropriate  fracture  mechanics  theory  for  concrete.  The  author  believes,  based  on  his 
experimental  results,  that  the  linear  elastic  fracture  mechanics  (LEFM)  model  is  still 
valid  for  the  modeling  of  crack  formation  and  propagation  in  concrete.  Although  my 
expertise  does  not  lie  in  this  particular  aspect  of  engineering,  I  tend  to  think  that  this 
formulation  will  produce  an  adequate  and  simple  model  for  the  cracking  and  fracture 
of  concrete  and  suggest  its  implementation  in  the  6.2  Nonlinear  Structural  Analysis 
Program. 

4.  Conclusions 

A  review  has  been  made  of  the  algorithms  developed  under  the  6.1  Structural 
Modeling  Program  to  evaluate  their  applicability  to  the  objectives  of  the  6.2  Nonlinear 
Structural  Analysis  Program.  The  major  conclusions  may  be  summarized  as  follows: 


The  proposed  reduced  coordinate  algorithms  should  have  a  solid  application  for 
the  problem  at  hand.  However,  it  is  recommended  that  the  "block  Lanczos"  method  be 
used  to  generate  the  approximating  subspace  of  Lanczos  vectors.  These  should  be 
obtained  taking  into  consideration  the  spatial  distribution  of  the  3  dimensional 
earthquake  components.  It  is  also  proposed  that  this  technique  be  used  in  the  context 
of  multi-level  substructuring.  Serious  consideration  should  be  given  to  the  fact  that  the 
programing  effort  involved  in  the  computer  implementation  of  the  proposed  Lanczos 
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algorithm  within  the  context  of  nonlinear  analysis  and  multilevel  substructuring  is 
considerable. 

The  algorithms  made  available  to  solve  the  nonlinear  contact  problem  are 
applicable  to  the  6.2  Program.  However,  it  is  necessary  that  simpler  contact  kinematic 
assumptions  be  made  for  the  3  dimensional  finite  element  implementation  of  these 
algorithms.  These  kinematic  hypotheses  could  be  improved  and  updated,  within  the 
finite  element  code,  to  more  complicated  schemes  as  research  in  the  area  of  the  3-D 
contact  problem  progresses. 

The  combination  of  the  boundary  element  and  the  finite  element  method  seems 
to  be  one  of  the  most  promising  approaches  for  the  solution  of  the  silent  boundary 
effect  in  the  time  domain.  Although  serious  consideration  is  being  given  to  it  in  current 
research,  no  general  algorithm  within  the  context  of  3  dimensional  elasto-dynamics  is 
ready  as  yet.  It  is  therefore  proposed  that  finite  element  techniques  in  the  time  domain, 
reported  in  the  literature,  be  used  in  conjunction  with  the  substructuring  procedure 
mentioned  above  to  account  for  the  energy  dissipation  at  the  boundaries  of  the  finite 
element  model. 

The  improved  bounding  surface  plasticity  model,  implemented  in  the  computer 
package  CLAY,  and  validated  through  experimental  testing,  should  be  of  direct 
application  to  model  the  nonlinear  conditions  of  the  soil  near  the  walls  of  the  drydock. 
The  linear  elastic  fracture  mechanics  (LEFM)  model  seems  to  be  valid  for  the  modeling 
of  crack  formation  and  propagation  in  concrete. 


Acknowledgment 

Support  for  this  work  has  been  provided  by  the  Naval  Civil  Engineering 
Laboratory  at  Port  Hueneme.  The  interest  of  Dr.  Shugar  and  Mr.  Holland  are  greatly 
appreciated. 


References 


1.  Mish  K.,  Romstad  K.,  Herrmann  L.,  "The  Solution  of  Nonlinear  Time  Dependent 
Problems  Using  Modal  Coordinates".  NCEL  Contract  Report  CR  86.004.  Port 
Hueneme.  California.  Dec  1985. 

2.  Mish  K.,  Herrmann  L.,  "The  Solution  of  Large  Time  Dependent  Problems  Using 
Reduced  Coordinates".  NCEL  Contract  Report  CR  87.010.  Port  Hueneme. 
California.  June  1987. 

3.  Landers  J.,  Taylor  R.  "An  Augmented  Lagrangian  Formulation  for  the  Finite 
Element  Solution  of  Contact  Problems".  NCEL  Contract  Report  CR  86.008.  Port 
Hueneme.  California.  March  1986. 

4.  Ju  J-W.,  Taylor  R.,  Cheng  L.  “A  Consistent  Finite  Element  Formulation  of 
Nonlinear  Frictional  Contact  Problems".  NCEL  Contract  Report  CR  87.008.  Port 
Hueneme.  California.  March  1987. 

5.  Cox  J.  V.,  Shugar  T.  A.  "A  Recursive  Boundary  Element  with  Initial  Error 
Investigation".  NCEL  Technical  Note  N-1730.  Port  Hueneme.  California.  Dec 
1985. 

6.  Herrmann  L.,  Kaliakin  V.,  Shen  C.  "Improved  Numerical  Implementation  of  the 
Bounding  Surface  Plasticity  Model  for  Cohesive  Soils".  NCEL  Contract  Report  CR 
86.004.  Port  Hueneme.  California.  Dec  1985. 

7.  Shen  CK.,  Zhu  ZY,  Herrmann  L.,  Kaliakin  V.,  "Validation  of  Bounding  Surface 
Plasticity  Theory  Using  Preliminary  Geotechnical  Centrifuge  Experiments".  NCEL 
Contract  Report  CR  86.011.  Port  Hueneme.  California.  July  1986. 

8.  Swartz,  S.  E.  "Applicability  of  Fracture  Mechanics  Methodology  to  Cracking  and 
Fracture  of  Concrete".  NCEL  Contract  Report  CR  86.006.  Port  Hueneme. 
California.  February  1986. 

9.  Bayo  E.,  "Numerical  Techniques  for  the  Evaluation  of  Soil-Structure  Interaction 
Effects  in  the  Time  Domain".  Ph.D.  Thesis.  University  of  California.  Berkeley. 
1983. 

10.  Bayo  E.,  Wilson  E.  "Finite  Element  and  Ritz  Vector  Techniques  for  the  Solution  to 
Three  Dimensional  Soil-Structure  Interaction  Problems  in  the  Time  Domain" 
Engineering  Computations,  Vol  1,  No  4,  pp  298-311.  Dec  1984. 

1 1 .  Lysmer  J.,  Udaka  T.,  Seed  H.B.  "FLUSH,  A  Computer  Program  for  Approximate 
3-D  Analysis  of  Soil  Structure  Interaction  Problems".  University  of  California 
Berkeley.  EERC  Report  75-30.  Nov  1975. 


12.  Lysmer  J.,  Raissi  M.,  Tajinian  F.,  Valdan  S.,  Ostadan  F.  "SASSI  -  A  System  for 
Analysis  of  Soil-Structure  Interaction".  University  of  California  Berkeley. 
Geotechnical  Report  #  81-02.  April  1981. 

13.  Day  S.  M.  "Finite  Element  Analysis  of  Seismic  Scattering  Problems".  Ph.  D. 
Thesis.  University  of  California.  San  Diego.  1977. 

14.  White  W.,  Valliapan  S.,  Lee  I.  K.  "Unified  Boundary  for  Finite  Dynamic  Models".  J. 
Eng.  Mech.  Div.  ASCE.  Vol  103,  pp  949-964.  1977. 

15.  Clough  R.,  Wilson  E.  "Dynamic  Analysis  of  Large  Structural  Systems  with  Local 
Nonlinearities".  Computer  Methods  in  Applied  Mechanics  and  Engineering.  Vol 
17-18,  pp  107-129.  1979. 

16.  Bathe  K.,  Gracewski  S."On  Nonlinear  Dynamic  Analysis  Using  Substructuring 
and  Mode  Superposition".  J.  Computer  and  Structures,  Vol.  13,  pp  699-707.  1981 

17.  Hurty  W.  C.,  Collins  J.  D.  "Dynamic  Analysis  of  Large  Structures  by  Modal 
Synthesis  Techniques".  Computer  and  Structures.  Vol  1,  pp  553-563,  1971. 

18.  Craig  R.R.  "Methods  of  Component  Mode  Synthesis".  The  Shock  and  Vibration 
Digest.  Vol  9,  pp3-10,  1977. 

19.  Guyan  R.J.  "Reduction  of  Stiffness  and  Mass  Matrices".  AIAA  Journal,  Vol  3  No  2. 
p  380.  1965. 

20.  Wilson  E.,  Bayo  E.  "Use  of  Special  Ritz  Vectors  in  Dynamic  Substructure 
Analysis".  ASCE  Journal  of  Structural  Engineering.  Vol  112,  No  8.  Aug  1986. 

21 .  Leger,  P.  "The  Use  of  Load  Dependent  Vectors  for  Dynamic  and  Earthquake 
Analyses".  Ph.  D.  Thesis.  University  of  California.  Berkeley.  1 986. 

22.  Nour-Amid  B.,  Clough  R.  "Block  Lanczos  Method  for  Dynamic  Analysis  of 
Structures".  Earth.  Eng.  Struct.  Dynam.  Vol  13,  pp  271-275.  1985. 

23.  Golub  G.H.,  Underwood  R.  "The  Block  Lanczos  for  Computing  Eigenvalues" 
Mathematical  Software  III,  edited  by  J.  Rice.  Academic  Press,  New  York,  pp  364- 
377.  1977. 

24.  Simon  H.,  "The  Lanczos  Algorithm  for  Solving  Linear  Systems".  Ph.  D.  Thesis, 
University  of  California.  Berkeley.  1982. 

25.  Nour-Amid,  B.,  Parlett  B.,  Taylor  R.,  "A  Newton-Lanczos  Method  for  Solution  of 
Non-Linear  Finite  Element  Equations".  Computers  and  Structures.  Vol  16,  No  1, 
pp  241-252.  (1983). 

26.  Nour-Amid  B.,  Clough  R.  "Dynamic  Analysis  of  Structures  Using  Lanczos 
Coordinates".  Earth.  Eng.  Struct.  Dynam.  Vol  12,  pp  565-577.  1984. 

14 


>*.  * 


.v.v yv . 


>■»«- L»  Itl  I.  I . 


27.  Ghaboussi  J.,  Wilson  E.  L.,  Isenberg  J.  "Finite  Element  for  Rock  Joints  and 
Interfaces".  Journal  of  Soil  Mechanics  and  Foundation  Division,  ASCE.  Vol  99. 
No  SM10.  pp833-848.  October  1973. 

28.  Hallquist,  J.O.  "A  Numerical  Procedure  for  Three  Dimensional  Impact  Problems." 
ASCE  Fall  Convention  and  Exhibit.  Preprint  2956.  October  1977. 

29.  Shugar  T.  A.,  Cox  J.  V.,  "An  Investigation  of  the  Indirect  Boundary  Element 
Method  in  One  and  Two  Dimensional  Elastostatics".  NCEL  Technical  Note  N- 
1664.  Port  Hueneme.  California.  May  1983. 


L 


WWW 


-V-l  - 


